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I.  INTRODUCTION 


1.1  BACKGROUND 

The  objective  of  the  Systems,  Science  and  Software  (S3) 
research  program  is  to  extend  our  present  understanding  of  the 
excitation  of  seismic  waves  by  underground  explosions  and 
earthquakes.  Toward  this  objective,  we  are  conducting  theo¬ 
retical  and  empirical  studies  of  ground  motions  from  the  two 
classes  of  sources.  In  particular,  our  efforts  are  directed 
toward  the  development  of  improved  methods  for  discriminating 
between  the  seismic  signals  from  earthquakes  and  explosions 
and  the  development  of  improved  methods  for  estimating  explo¬ 
sion  yield. 

This  report  summarizes  the  work  done  during  the  second 
three-month  period  of  the  contract. 


1.2  SUMMARY  OP  RESEARCH  DURING  THIS  QUARTER 

Our  work  during  this  quarter  has  included  research  in  a 
number  of  areas.  Research  projects  are  briefly  summarized 
below. 

Source  Calculations 

A.  Calculations  of  Nuclear  Explosions  in  Granite 

One-dimensional  and  two-dimensional  finite  difference 
calculations  of  the  PILEDRIVER  event  have  been  made  and  are 
in  excellent  agreement  with  the  observations.  A  brief  dis¬ 
cussion  of  these  calculations  is  given  in  Section  V.  Also 
discussed  in  that  section  are  some  one-dimensional  calcula¬ 
tions  of  the  French  explosions  in  the  Sahara.  In  Section  VI 
we  discuss  some  recent  calculations  of  decoupled  explosions 
in  salt.  The  results  are  in  good  agreement  with  the  SALMON/ 
STERLING  experiment. 
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B.  Earthquake  Modeling  on  the  ILLIAC  IV  Computer 

Our  earthquake  calculations  done  on  the  ILLIAC  computer 
have  so  far  included  only  a  rather  simple  model  for  the  physics 
of  earthquake  faulting.  The  great  advantage  of  the  TRES  code 
on  the  ILLIAC  computer  is  that  realistic  physics  can  be  incor¬ 
porated  into  the  model.  In  Section  IV  we  describe  a  recently 
developed  algorithm  for  a  three-dimensional  numerical  treat¬ 
ment  of  fracture  on  a  fault  plane.  This  model  is  now  being 
incorporated  into  the  TRES  program  on  the  ILLIAC. 

Discrimination 

During  this  reporting  period  we  completed  the  MARS  pro¬ 
cessing  for  127  events  recorded  at  the  Priority  1  stations. 

After  several  attempts,  and  a  recopy  of  the  original  tape,  we 
have  been  unable  to  convert  data  for  event  149.  In  addition, 
we  have  not  received  digital  seismograms  for  events  274  through 
278  for  any  of  the  Priority  2  stations. 

Preliminary  variable  frequency  magnitude  (VFM)  results 
for  the  Priority  1  stations  were  presented  at  the  VSC  meeting 
on  April  11,  1979  in  Alexandria,  Virginia.  The  degree  of 
separation  of  earthquakes  and  presumed  explosions  appears  to 
be  dependent  on  several  factors.  In  addition  to  possible 
spectral  differences  in  the  earthquake  and  explosion  source 
functions,  the  propagation  path  in  the  vicinity  of  a  particular 
station  and  the  level  of  background  earth  noise  both  contribute 
to  this  separation.  To  date  the  station  RKON,  sited  on  the 
Canadian  shield,  has  yielded  the  maximum  separation  of  the  pre¬ 
sumed  explosions  from  the  majority  of  earthquakes.  In  contrast, 
the  two  event  populations  exhibit  the  least  amount  of  separation 
at  the  station  ANMO,  located  within  the  Basin  and  Range  Province. 
Two  reports,  describing  detailed  results  from  the  Priorities  1 
and  2  station  sets,  are  in  preparation. 
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Small-Scale  Experiments 

Recent  progress  in  simulating  underground  explosions  by 
using  small  PETN  charges  in  concrete  blocks  is  described  in 
Section  II. 

Research  Activities  at  the  Reston  Geophysics  Office 

The  research  project  includes  several  tasks  being  done 
by  J.  R.  Murphy  and  his  colleagues  in  our  Reston  Office.  Pro¬ 
gress  on  these  tasks  is  summarized  in  Section  III. 
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II.  SEISMIC  MODELING  EXPERIMENTS 


2.1  INTRODUCTION 

The  objective  of  these  experiments  will  be  to  make 
measurements  of  free  surface  motion  due  to  the  detonation  of 
one  or  more  small  high  explosive  (HE)  charges  within  a  uniform, 
well-characterized  medium.  The  data  will  be  used  as  bench¬ 
marks  for  the  finite  difference  codes  used  to  predict  the 
seismic  motions  generated  by  underground  explosions.  The 
principal  activity  to  date  has  been  the  development  and  test¬ 
ing  of  the  HE  charges;  see  the  discussion  below.  The  initial 
experiment,  designed  to  check  instrument  performance  and  con¬ 
firm  previous  measurements,  will  be  conducted  by  13  April  1979. 


2.2  HE  CHARGE  DEVELOPMENT 

There  are  four  main  requirements  for  the  charges: 

1.  They  should  be  spherical. 

2.  The  power  and  energy  output  should  be 
uniform  from  charge  to  charge. 

3.  The  performance  of  the  explosive  should 
be  well  characterized  so  that  the  codes 
cam  accurately  model  the  detonation. 

4.  The  charges  should  be  small  so  that  a 
reasonably  sized  block  of  the  medium 
(of  order  one  cubic  meter)  gives  elastic 
motions  on  its  free  surfaces. 

A  few  tenths  gram  of  PETN,  initiated  by  am  exploding 
bridgewire  and  contained  within  a  plastic  shell  of  one  centi¬ 
meter  diameter,  meets  all  these  requirements.  Figure  1  shows 
some  details  of  construction.  The  shell  is  formed  from  two 
hemispheres.  The  bridgewire  and  its  leads  are  sealed  into 
one  hemisphere.  A  small  loading  port  is  cut  from  the  other 
hemisphere.  After  the  hemispheres  are  glued  together,  a 
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Figure  1.  Spherical  explosive  charge  (section  view,  10  x  scale) 
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measured  amount  of  PETN  powder  is  poured  into  the  shell  and 
then  the  port  is  sealed  in  place.  All  seals  are  made  with 
styrene  solvent  cement.  This  assembly  technique  ensures 
accurate  centering  of  the  bridgewire  and  good  waterproofing 
for  the  explosive. 

We  evaluated  the  performance  of  these  charges  by  deto¬ 
nating  them  in  a  small  tank  of  water  and  photographing  the 
tests  with  a  framing  camera  operating  at  about  500,000  frames 
per  second.  In  order  to  see  both  the  expansion  of  the  HE  and 
the  shock  wave  in  the  water,  the  tank  was  backlit.  Such  photo¬ 
graphs  are  essentially  shadowgraphs*,  thus  the  ideal  light 
source  is  a  distant,  bright  point  emitter.  We  approximated 
the  ideal  with  a  commercial  photo  flash  "strobe"**  located 
about  one  tenth  meter  behind  the  four  inch  cubical  water  tank. 
The  film  was  a  high  speed,  high  contrast  type  (Kodak  2475,  ASA 
1,000)  and  the  effective  aperature  of  the  camera  was  f/28. 

We  also  placed  a  pressure  sensor  20.3  mm  from  the 
surface  of  the  HE  sphere  to  measure  the  water  shock.  The 
sensor  was  a  bar  gauge  type.  The  actual  probe  in  the  water, 
a  one-eight  inch  diameter,  100  mm  long  piece  of  oil-hardened 
steel  drill  rod,  acoustically  transmits  the  pressure  on  the 
end  of  the  rod  to  the  transducer,  x-cut  quartz. 

Some  examples  of  the  photographs  for  the  first  test  on 
14  March  1979  are  shown  in  Figure  2.  Because  the  camera  cannot 
by  easily  synchronized  with  the  high  voltage  detonating  pulse, 
absolute  time  is  uncertain  for  each  frame.  Defining  frame  1 
as  showing  first  evidence  of  detonation  leads  to  the  timing 
shown  in  the  figure.  Frame  -1  shows  the  unexploded  sphere 


4* 

The  mean  size  of  the  PETN  particles  must  be  small  compared  to 
the  diameter  of  the  bridgewire  (25  um)  to  ensure  high  order 
detonation  of  the  explosive. 

*The  HE  products  are  at  a  temperature  about  2500  to  3000°K  which 
makes  them  weakly  self  luminous  in  the  red  and  infrared  parts 
of  the  spectrum. 

**Rated  guide  number  of  60  for  a  film  speed  of  ASA25 . 
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Frame  2 
t*3 . 7  js 


Frame  7 
t*13 . 1-s 


•e  2. 


Framing  camera  photographs  of  detonation  o 
sphere  in  water,  first  test  of  14  March  19 
Arrow  shows  shock  front.  (Scale  x  0.88) 
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(painted  black)  supported  on  its  lead  wires  in  front  of  a  10  ram 
spacing  grid  pattern.  Expansion  of  HE  gases  is  evident  even  in 
frame  1.  Measurements  from  these  prints  give  an  initial  water 
shock  velocity  of  2.8  mm/ys  and  a  maximum  particle  velocity 
(i.e.,  the  HE  gas  motion)  of  0.6  to  0.7  ram/us.  Both  velocities 
are  consistent  with  an  initial  water  shock  pressure  of  1.8  GPa 
(18  kbar) .  Based  on  a  packing  density  of  0.8  g/cc,  the  equation 
of  state  for  PETN  (N.  Rimer,  private  communication)  also  pre¬ 
dicts  an  initial  pressure  of  1.8  GPa  in  the  HE  products.  In 
frame  3,  the  thin  brigh  band  along  the  sphere's  equator  is 
presumably  due  to  the  hot  HE  products;  the  joint  between  the 
hemispheres  is  probably  thinner  than  the  rest  of  the  shell. 

In  frame  7  the  shadow  at  the  top  shows  that  the  water  shock 
has  moved  a  few  millimeters  beyond  the  end  of  the  bar  gauge. 

By  this  time,  the  shock  has  slowed  to  1.7  mm/ys*  and  the  HE 
products  are  expanding  at  0.4  mm/ys. 

Defining  asymmetry  in  percent  as 

[Polar  Diameter]  -  [Equatorial  Diameter]  ..  ,  „„ 
- tPolar"  Diameter]  -  x  100 

the  HE  products  show  a  maximum  of  three  percent  asymmetry  and 
the  shock,  minus  five  percent,  for  this  test.  By  10  ys,  both 
measures  were  less  than  one  percent.  Other  tests  show  similar 
or  smaller  asymmetries.  The  centers  of  the  sphere,  expanding 
HE  gases,  and  water  shock  are  coincident  to  within  a  millimeter. 

Figure  3  gives  a  typical  pressure  record.  The  finite 

risetime  (1  ys)  and  ringing  after  the  peak  (period  of  about 

+ 

2  ys)  are  characteristics  of  a  bar  subjected  to  a  step  func¬ 
tion  pressure  input.  This  record  implies  a  true  shock  pressure 

* 

Sonic  velocity  in  water  is  1.65  mm/ys. 

^Dispersion  in  the  propagation  of  compressional  waves  along 
the  bar  account  for  these  features. 
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of  0.2  GPa  (2  kilobar)  with  a  submicrosecond  rise  time  and  1/e 
decay  of  about  5  us.  These  data  are  consistent  with  the  veloc¬ 
ity  data  from  the  photographs. 

As  a  result  of  these  tests,  we  can  now  generate  con¬ 
sistent,  symmetric,  high-order  detonations  in  small  explosive 
charges.  In  anticipation  of  our  planned  seismic  modeling 
experiments,  we  will  also  photograph  the  simultaneous  explosion 
of  two  or  three  charges  in  a  water  tank  to  confirm  equal  detona¬ 
tion  rates  between  charges. 
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III.  SUMMARY  OF  RESEARCH  ACTIVITIES  AT  THE 
RESTON  GEOPHYSICS  OFFICE 

3.1  MAGNITUDE- YIELD  IMPROVEMENT 

An  analysis  effort  has  been  initiated  to  assess  the 
detectability  of  decoupled  explosions  at  regional  distances 
as  a  function  of  source  conditions,  epicentral  distance, 
frequency  and  background  noise  conditions.  The  investigation 
is  focusing  on  a  critical  evaluation  of  the  model  described  by 
Rodean  of  LLL  at  the  recent  ARPA  Decoupling  Conference.  An 
attempt  is  being  made  to  place  reasonable  upper  and  lower 
bounds  on  each  of  the  elements  which  appear  in  this  model  in 
order  to  provide  estimates  on  the  range  of  the  expected  fre¬ 
quency  dependent  detectability  at  regional  distances.  For 
example,  the  source  model  used  in  Rodean 's  simulation  is 
essentially  the  original  Latter  theoretical,  low  frequency 
approximation  corresponding  to  a  simple  step  in  pressure  acting 
on  the  cavity  wall.  Experience  has  shown  that  this  provides  a 
lower  bound  estimate  on  the  amplitude  of  the  decoupled  source 
function  (i.e.,  it  predicts  a  decoupling  factor  of  200  for 
Sterling  versus  the  inferred  empirical  value  of  about  70) . 

The  corresponding  upper  bound  estimate  for  the  source  function 
has  been  obtained  by  reducing  the  Latter  low  frequency  decou¬ 
pling  to  match  Sterling  experience  and  modifying  the  high  fre¬ 
quency  spectral  composition  to  include  the  effect  of  the  initial 
pressure  spike  on  the  cavity  wall.  Similar  lower  and  upper 
bounds  are  being  defined  for  the  controlling  variables  of  the 
propagation  path  and  local  site  noise  conditions  for  use  in  a 
series  of  bounding  detectability  simulations. 


3.2  GROUND  MOTION  ANALYSIS 

Effort  on  this  task  has  focused  on  the  analysis  of  the 
free-field  data  from  the  Merlin  event.  The  data  recorded  at 
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shot  depth  from  this  event  over  the  distance  range  from  about 
200  to  750  m  are  unusual  because  of  the  presence  of  a  late¬ 
time  (i.e.,  greater  than  1.2  seconds),  relatively  long-period 
(-o  0.8  second)  secondary  arrival.  While  this  waveform  is  not 
particularly  noticeable  on  accelerograms  it  is  detectable  on 
the  velocity  records  and  becomes  quite  prominent  on  the  dis¬ 
placement  records  where  amplitudes  approach  those  of  the  direct 
arrivals.  Three  possible  sources  of  this  arrival  have  been 
considered:  (1)  free  surface  reflected  phases  (i.e.,  pP  or  pS) , 

(2)  reflections  from  a  deep  interface,  and  (3)  spall  closure. 

First,  with  regard  to  surface  reflected  phases,  the 
expected  arrival  times  for  pP  (0.5  to  0.8  seconds)  and  pS  (0.7 
to  1.0  seconds)  are  significantly  earlier  than  the  observed 
arrival  time  of  1.3  to  1.7  seconds.  Moreover,  the  observed 
radial  displacement  amplitudes  associated  with  this  arrival 
are  nearly  constant  over  the  range  200  to  750  m,  which  is 
inconsistent  with  the  behavior  expected  for  pP  or  pS. 

The  second  possibility  considered  was  that  the  arrival 
represents  a  reflection  from  a  deep  interface.  However,  in 
order  to  explain  the  late  arrival  time,  the  interface  would 
have  to  be  located  about  1  km  below  the  shot  and  the  apparent 
velocity  across  the  shot  depth  stations  from  a  reflection  at 
this  depth  would  be  expected  to  be  much  greater  than  the 
observed  apparent  velocity  of  about  1800  m/sec.  Moreover,  the 
amplitude  of  a  reflection  from  this  depth  would  not  be  expected 
to  be  very  large  on  the  radial  component. 

Having  eliminated  surface  reflections  and  deep  reflec¬ 
tions  as  explanations  for  the  observation,  we  are  left  with 
spall  closure.  It  is  noted  that  the  surface  ground  motion 
records  from  Merlin  also  show  a  relatively  large  arrival  at 
late  time  (0.9  to  1.3  seconds)  for  stations  with  slant  ranges 
from  about  300  to  700  m.  These  signals  coincide  approximately 
with  the  arrival  of  what  appears  to  be  a  rather  prolonged  spall 
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closure  signal,  although  the  observed  period  is  substantially 
longer  than  that  of  the  impulsive  motion  generally  observed  for 
spall  closure  at  near  stations.  Calculations  indicate  that  a 
P  wave  generated  near  the  surface  at  the  times  and  positions 
indicated  by  the  surface  recordings  would  arrive  at  the  source- 
depth  stations  at  about  the  same  time  as  the  unusual  observed 
signals. 

In  fact,  quantitative  simulation  of  the  spall  closure 
as  an  axisymmetric  distribution  of  vertical  forces  on  the  sur¬ 
face  produces  source  depth  pulses  which  agree  remarkably  well 
with  the  observations  with  regard  to  arrival  time  and  amplitude/ 
distance  dependence.  However,  the  sense  of  the  predicted 
motion  is  out  with  respect  to  the  explosive  source  as  opposed 
to  the  observed  inward  displacement  associated  with  the  pulse 
in  question.  Near-field  analytical  solutions  are  currently 
being  examined  in  an  attempt  to  explain  this  puzzling  phase 
problem. 
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IV.  A  FRACTURE  CRITERION  FOR  FINITE  DIFFERENCE 
MODELING  OF  EARTHQUAKE  FAULTING 

4.1  INTRODUCTION 

In  our  previous  three-dimensional  simulations  of  earth¬ 
quake  faulting,  we  have  modeled  faulting  as  a  propagating 
stress  relaxation  on  a  plane.  We  have  bypassed  much  of  the 
physics  of  the  rupture  process,  however.  Specifically,  we 
have  (1)  prescribed  the  propagation  of  the  fault  surface,  (2) 
constrained  slip  to  be  in  a  prescribed  direction  (parallel  to 
the  equilibrium  shear  stress  in  the  fault  plane) ,  and  (3)  re¬ 
quired  the  fault  to  heal  permanently  once  the  slip  reversed 
direction. 

We  have  developed  a  new  algorithm  governing  three- 
dimensional  frictional  sliding  which  permits  a  more  physical 
treatment  of  the  fault,  relaxing  constraints  2  and  3.  Changes 
in  direction  of  sliding,  as  well  as  arrest  and  restarting  of 
slip,  are  permitted  during  the  faulting  process.  Arrest  of 
sliding  is  determined,  in  the  new  algorithm,  by  fault  dynamics, 
rather  than  being  somewhat  arbitrarily  enforced  when  the  slid¬ 
ing  direction  reverses.  In  fact,  the  notion  of  a  velocity 
reversal,  which  is  ill-defined  when  the  slip  direction  is  not 
specified  a  priori,  is  dispensed  within  the  new  algorithm. 

The  new  frictional  sliding  algorithm  is  ideally  suited 
for  incorporating  a  fracture  criterion,  so  that  rupture  ad¬ 
vance  is  governed  by  rock  strength  rather  than  being  arbitrarily 
prescribed  (so  that  we  can  relax  constraint  1) .  A  physically 
acceptable  fracture  criterion  should  (a)  dissipate  energy  at 
the  crack  tip  and  (b)  lead  to  bounded  stresses  and  velocities. 

A  third  requirement  is  motivated  by  numerical  considerations; 
that  is,  an  acceptable  fracture  criterion  should  be  able  to  be 
formulated  so  as  to  be  independent  of  zone-size,  or  nearly  so. 
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The  first  two  requirements  are  satisfied  by  the  slip¬ 
weakening  model  investigated  by  Ida  (1972)  and  Andrews  (1976a) . 
A  two-dimensional  numerical  study  by  Andrews  suggests  that  it 
is  also  zone-size  independent.  This  model  of  failure  presumes 
that  the  frictional  traction,  at  a  given  point  on  the  prospec¬ 
tive  fault  plane,  is  a  prescribed  function  of  the  total  slip 
which  has  occurred  at  the  point.  This  can  be  considered  an 
idealization  of  a  strain-weakening  law. 

In  the  next  section  we  will  outline  the  three-dimensional 
treatment  of  frictional  sliding  in  the  context  of  a  simple 
Coulomb  law.  Then,  the  incorporation  of  a  slip-weakening  law 
requires  only  a  minor  refinement  of  the  Coulomb  law  and  is 
described  in  the  subsequent  section. 


4.2  FINITE  DIFFERENCE  TREATMENT  OF  COULOMB  FRICTION  IN 

THREE-DIMENSIONAL  CODES 

A  difficulty  arises  in  three-dimensional  numerical 
schemes  when  treating  Coulomb  frictional  sliding  on  an  inter¬ 
ior  surface,  due  to  the  fact  that  the  sliding  direction  is, 
in  general,  not  known  a  priori.  A  solution  is  presented  here 
which  resolves  the  difficulty.  The  method  has  been  success¬ 
fully  incorporated  into  the  TRES  three-dimensional  finite 
difference  code. 

Physically,  the  frictional  traction  on  a  slip  surface 
is  required  to  contribute  an  acceleration  in  the  direction 
opposed  to  the  instantaneous  slip  velocity.  However,  in 
numerical  time-stepping  schemes,  velocity  and  acceleration 
are  centered  differently  in  time.  This  is  not  a  trivial 
discrepancy;  if  the  direction  of  the  frictional  acceleration 
contribution  at  time  t  is  determined  from  the  velocity  at 
time  t  -  At/2  (where  At  is  the  time  increment) ,  friction  will 
generally  have  a  component  in  the  direction  of  sliding.  Even 
if  this  component  is  initially  small,  it  will  add  energy  to 
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the  system,  and  the  resulting  behavior  is  unstable.  A  re¬ 
lated  difficulty  is  that  since  the  slip  has  two  components, 
it  becomes  impossible  to  detect  a  slip  velocity  zero,  since 
velocity  is  known  only  at  discrete  times.  Thus,  it  is  not 
obvious  how  to  formulate  a  criterion  for  the  arrest  of  slip 
at  a  point. 

We  assume  that  slip  is  confined  to  the  plane  z  »  0, 
and  is  governed  by  a  Coulomb  law.  £(t)  denotes  the  slip  vec¬ 
tor  (+z  side  relative  to  -z  side)  at  time  t  at  a  point  in  the 
plane  z  =  0.  Only  tangential  slip  is  permitted,  t_  denotes 
the  tangential  traction  exerted  on  the  +z  half space  at  z  *  0. 
t  is  the  value  of  this  tangential  traction  required  to  en- 
force  continuity  of  velocity,  and  is  the  magnitude  of  the 
sliding  frictional  stress  (af  is  positive  and  proportional 
to  the  compressive  stress,  -  a  ,  on  the  slip  plane).  With 
these  definitions,  the  friction  law  takes  the  form  of  the 
following  boundary  condition  on  z  »  0: 


if  llj  *  of 
if  llj  <  of 


(1) 


We  write  the  x  and  y  components  of  displacement  u  at 
the  slip  plane  z  *  0,  as  the  sum  of  a  continuous  part  (q) 
and  a  discontinuity  (s) : 


u 


q  +  s 
Mx  x 


u. 


q  +  s 

My  y 


16 


■S.I 


SYSTEMS.  SCIENCE  AND  SOFTWARE 


(u_  is  continuous).  Assuming  an  explicit  difference  scheme, 
the  $'s  and  s's  are  determined  at  times  (n-l/2)At  (At  is  the 
time  step,  n  is  an  integer)  and  the  s's  and  q's  at  times 
nAt  by  recursion  relations  of  the  form 

•  n  .n-1  ...  f  _  /n-1  n-1  n-1  n-1  n-l\  ,  _ 

S  38  s  +  At  F  ( s  ,s  ,g  ,q  ,u  l+T 

X  x  [  x\  X  y  4x  4y  '  uz  )  > 

•  n  .n-1  ,  . .  [  _  /..n-1  n-1  n-1  n-1  n-1  \  ,  _ 

sy  y  +  4t  [  Fy\sx  -  ay  •  ’’x  -  «y  -  uz  )  +  Ts 

4"  -  q"-1  .its  (s"'1,  sn'X,  q"'1,  q""1,  un‘l) 

4n  =  q11-1  +  At  G  Is n"1,  s*"1,  q11"1,  q*"1 ,  u^1) 

Hy  Hy  y\  X  '  y  '  Hx  '  Hy  ’  uz  / 

•n  _  .n-1  _  /n-1  n-1  n-1  n-1  n-l\ 

u_  =  u_  +  At  Gls  ,  s„  ,  q„  ,  q„  ,  u_  ) 

zz  zyx  y  x  y  z  / 


n  n-1  ...  *n 
3x  *  3x  +  4t  sx 


•5*1  - 


n-1  ,  . .  •  n 

qx  +  At  qx 


n-1  .  ..  • n 
+  At  qy 


»rl  +  **  a; 


F  ,  F  ,  G  ,  G  ,  and  G  are  spatial  difference  operators,  the  G’s 
x  y  x  y  z 

representing  continuous  accelerations,  the  F's  a  discontinuity 
in  acceleration. 


svsrrMS.  scun cc  and  software 


T  and  T  represent  the  acceleration  due  to  the  boundary 
x  y 

stresses  from  Equation  (1) .  The  most  obvious  form  for  Tx  and 
Ty  would  be 


(A  is  positive  constant,  involving  density  and  spatial  zone 
size,  which  converts  stress  to  nodal  acceleration) . 

Unfortunately,  Equation  (3) ,  combined  with  Equations  (2a) 
and  (2b) ,  leads  to  an  instability,  due  to  the  fact  that  the  T's 
are  centered  (1/2) At  earlier  in  time  than  the  F's.  A  solution 
proposed  by  Day  (1977)  has  been  found  to  provide  a  stable  solu¬ 
tion.  Centering  is  modified  by  replacing  Equation  (3)  with  an 
average  over  two  time  steps: 
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However,  a  new  definition  of  Tc  associated  with  Equation  (5) 
is  necessary,  and  will  be  derived  below. 

i  Inserting  Equation  (5)  into  2a  and  2b  leads  to  a 

pair  of  coupled  nonlinear  equations  for  unknowns  s"  and  s”: 
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MMsJf 
N2  *  (s)f 


where 


sn_1  +  AtF 

X  X 


AtAa^ 


•  n- •  1 

3y  +  AtFy 


AtA  o, 


Vr  •  (r)f 

_ s^1 


AtACJ, 


These  have  the  solution  (see  Day,  1977) 


•n  _ 
sx  *  Cx 


sn  =  C 

y  y 


2  2 ' 

C  +  C 
x  y , 


2  2 

C  +  C 
x  y 


if  (cf 


if  C  +C  2  B 
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which  provides  a  stable  scheme  for  updating  s  and  s  .  If 

(2  2\^/2  x  y 

C1  +  C2/  <  B'  no  solution  exists  for  Equation  (6),  which 

suggests  that  we  redefine  Tc  as 


and  replace  Equation  (2a) ,  (2b) ,  (3a) ,  and  (3b)  with 


Equation  (8)  provides  well-behaved  simulation  of  Coulomb 
friction.  The  condition  Tc  <  o^,  with  Tc  given  by  Equation 
(7),  is  the  criterion  for  arrest  of  sliding.  Mathematically, 
this  condition  says  that  the  velocity  must  be  continuous  when 
no  solution  exists  to  the  system  (6) .  Physically,  it  says 
that  slip  must  be  arrested  at  time  nAt  if  a  slip  velocity 
zero  will  occur  due  to  frictional  deceleration  during  the 


subsequent  (1/2) At.  If  the  slip  velocity  were  not  set  to  zero, 
the  frictional  force  would  be  acting  to  accelerate  slip  during 
part  of  that  (1/2) At. 

The  TRES  finite  difference  code  (Cherry,  1977)  has  pre¬ 
viously  avoided  the  problem  by  artificially  constraining  slip 
to  occur  in  a  prescribed  direction  only.  Equation  (8)  has 
been  incorporated  into  TRES,  permitting  this  constraint  to  be 
relaxed.  The  algorithm  governing  frictional  sliding  in  TRES 
has  actually  been  considerably  simplified  as  a  result. 
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4.3 


SLIP-WEAKENING  CONSTITUTIVE  MODEL 


Ida  (1972)  investigated  the  two-dimensional,  anti¬ 
plane-strain  dynamics  of  a  slip-weakening  failure  model.  In 
this  model,  a  finite  shear  strength,  aQ(0),  is  assigned  to  the 
fault  prior  to  initiation  of  sliding.  Slip  commences  at  a 
point  when  necessary  to  prevent  a  stress  concentration  in 
excess  of  Cq( 0) .  The  relative  displacement,  with  amplitude 
denoted  by  s,  is  assumed  to  weaken  the  fault  plane  at  that 
point  according  to  some  prescribed  function  ctq(s)  until  the 
fault  loses  all  cohesion  and  reaches  the  kinetic  friction 
level,  aQ( D) ,  where  D  is  the  slip  required  to  eliminate 
cohesion.  The  amount  of  energy  required  per  unit  area  to 
destroy  cohesion,  denoted  2G,  is 


D 

-/  - 

•  A 


Q(s)  -  cr0(D)  ]  ds 


A  simple  slip-weakening  law  can  be  immediately  intro¬ 
duced  in  Equation  (8) .  All  that  is  necessary  is  that  af  in 
Equation  (8)  be  made  a  function  of  the  slip: 


Of  = 


a0(s). 


VD>' 


if  s  <  D 


if  s  >  0 


where  D  is  a  specified  constant,  and  the  function  aQ(s)  is 
a  prescribed  function  for  0  <_  s  £  D.  Ida's  analysis  suggests 
that  if  Cq  is  assumed  to  be  fairly  smooth,  the  behavior  of 
the  slip  function  is  insensitive  to  the  precise  shape  of  cQ. 
Currently,  we  prescribe  to  be  linear. 
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An  algorithm  based  on  Equation  (8) ,  with  determined 
as  in  Equation  (10) ,  bus  been  communicated  to  the  Institute 
for  Advanced  Computation  (IAC)  programmers.  They  are  in  the 
process  of  incorporating  the  model  into  I4TRES ,  the  ILLIAC  IV 
version  of  TRES.  When  the  ILLIAC  program  modifications  are 
completed,  we  will  begin  examining  the  role  of  the  surface 
energy,  G,  in  determining  the  shape  and  rate  of  advance  of  the 
rupture.  Of  particular  interest  initially  will  be  whether 
rupture  velocities  exceeding  the  S  wave  speed  develop  for  three- 
dimensional  rupture  models,  as  they  persistently  do  for  two- 
dimensional  numerical  models  (Andrews,  1976b?  Das  and  Aki, 

1977) .  Subsequently,  we  will  look  at  the  role  of  stress 
heterogeneities  in  stopping  the  rupture. 

This  constitutive  model  can  be  viewed  as  an  idealiza¬ 
tion  of  a  strain-weakening  constitutive  law.  Andrews  (1976a) 
has  argued  that  G  cannot  be  considered  a  material  constant, 
at  least  once  the  fault  has  grown  beyond  a  certain  limiting 
size.  Eventually,  inelastic  dissipation  is  expected  to  in¬ 
volve  a  widening  zone  about  the  fault  plane.  At  this  point 
the  slip-weakening  idealization  becomes  inadequate,  unless  G  is 
increased  to  properly  account  for  the  increasing  inelastic 
work.  Our  approach  is  to  analyze  the  simple  slip-weakening 
model,  recognizing  its  probable  limitations,  and  then  general¬ 
ize  it  to  include  realistic  inelastic  behavior  off  the  fault 
plane . 
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V.  SOURCE  CALCULATIONS 


An  extensive  one-  and  two-dimensional  effort  is  con¬ 
tinuing  to  simulate  the  ground  motion  from  the  PILEDRIVER  event 
in  NTS  fractured  granodiorite.  Of  particular  interest  are  the 
effects  of  the  free  surface  (spall,  overburden,  etc.)  on  the 
far-field  ground  motion.  The  two-dimensional  calculations  to 
be  believable  should  be  in  agreement  with  the  measured  free- 
field  and  free  surface  velocity  data  (Perret,  1968;  Hoffman 
and  Sauer,  1969) .  It  is  well  known  that  these  data  cannot  be 
matched  using  the  triaxial  failure  envelope  measured  in  the 
laboratory  either  for  competent  or  fractured  granite.  We  have 
reported  earlier  (Bache,  et  al . ,  1978)  the  results  of  one¬ 
dimensional  free  field  calculations  made  using  an  effective 
stress  model  which  accounts  for  the  water  present  in  the  jointed 
granite.  These  calculations  were  in  excellent  agreement  with 
the  measured  cavity  radius  and  in  fairly  good  agreement  with 
the  velocity  data  from  two  shot  level  stations  reported  by 
Perret.  However,  a  preliminary  two-dimensional  calculation 
made  with  the  same  material  properties  did  not  give  a  good 
match  to  the  rest  of  the  free-field  data  or  to  the  free  sur¬ 
face  data.  A  longer  duration  of  the  positive  velocity  pulse 
was  required  to  better  match  the  free  field  data. 

A  series  of  one-dimensional  free  field  calculations 
has  been  made  to  investigate  the  effect  of  changes  at  the 
low  pressure  end  of  the  failure  surface  (the  higher  pressure 
failure  envelope  for  fractured  granite  is  well  defined) .  For 
a  given  choice  of  Yq,  the  strength  at  P  *  0  (see  Bache,  et  al. , 
1978  for  details  of  the  shear  failure  model) ,  we  attempted  to 
match  the  calculation  to  the  free  field  data  through  a  judi¬ 
cious  choice  of  Pc,  the  pressure  at  which  all  air-filled  pores 
are  closed.  Our  strength  model  directly  links  the  air-filled 
porosity  to  the  effective  pressure,  the  difference  between  the 
pressure  in  the  rock  and  the  pressure  in  the  pore  water.  When 
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p  is  reached,  the  rock  is  considered  to  be  in  pressure  equi- 
c 

librium  so  that  the  effective  pressure  is  zero.  Since  the 
shear  strength  of  the  rock  is  a  function  of  the  effective 
stress,  a  high  value  of  Pc  means  that  a  smaller  region  sees 
the  lower  strength.  We  believe  that  Pc  cannot  be  more  than 
a  few  kilobars  as  most  of  the  porosity  for  the  granodiorite 
is  considered  to  be  in  the  form  of  easily  closed  cracks. 

The  velocity  gauge  data  at  the  two  Perret  shot  level 
stations  were  successfully  matched  in  two  free  field  calcula¬ 
tions.  Figures  4  and  5  show  the  agreement  between  the  data 
and  Calculation  1  made  with  a  relatively  low  strength  YQ  = 

20  bars  and  with  P  =  2.0  kb.  Figures  6  and  7  show  compari- 

c 

sons  between  the  data  and  Calculation  2  for  Y  =  100  bars  and 

o 

Pc  =  1.0  kb.  There  appears  to  be  little  to  choose  between  the 
two  calculations  based  on  the  velocity  data  or  even  based  on 
the  transform  of  the  reduced  velocity  potential  (RVP) .  (Cal¬ 
culation  1  has  a  slightly  lower  corner  frequency  but  a  larger 
static  RDP . )  However,  the  different  values  of  Yq  should 
strongly  influence  the  yielding  and  spall  near  the  free  sur¬ 
face.  A  choice  may  be  made  based  on  the  calculated  cavity 
radii;  47.4  meters  for  Calculation  1  versus  39.5  meters  for 
Calculation  2.  The  more  widely  accepted  measurements  of  cavity 
radius  are  40  meters,  based  on  drillback  into  the  lower  half  of 
the  PILEDRIVER  cavity,  and  44.5  meters  calculated  from  the  void 
volume  of  the  chimney  determined  from  chimney  pressurization 
tests.  (The  latter  number  is  in  better  agreement  with  the 
scaled  cavity  radius  from  HARDHAT,  a  nuclear  event  in  the  same 
material  as  PILEDRIVER. )  A  final  determination  on  which  cal¬ 
culation  to  accept  has  not  been  made  at  this  time  and  will,  of 
course,  depend  upon  the  results  of  two-dimensional  studies. 

Several  additional  coarsely  zoned  two-dimensional  cal¬ 
culations  have  been  made,  each  one  incorporating  more  of  the 
geologic  structure  of  the  PILEDRIVER  site.  Calculation  PD15 
used  the  time  of  arrival  data  as  reported  by  Hoffman  and  Sauer 
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PARTICLE  VELOCITY  (10  ft/sec) 


TIME  (10*1  sec) 


Figure  4.  Comparison  between  data  and  Calculation  1 

(Pc  *  2.0  kb,  Y0  =  20  bars)  at  first  Perret 
station  (668  feet  from  WP) . 
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TIME  (10  sec) 


Figure  5 


Comparison  between  data  and  Calculation  1 
(Pc  *  2.0  kb,  Y0  *  20  bars)  at  second  Perret 
station  (1543  feet  from  WP) . 
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TIME  UO”1  sec) 


Figure  6.  Comparison  between  data  and  Calculation  2 

(Pc  ■  1.0  kb,  Y0  »  100  bars)  at  first  Perret 
station. 
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to  locate  two  horizontal  layers  approximately  350  and  420  meters 
above  the  working  point.  This  calculation  used  the  material 
properties  of  Calculation  1  for  the  working  point  layer  of 
material. 

Another  two-dimensional  calculation  (PD16)  used  the 
material  properties  of  Calculation  2  for  the  working  point 
material  and  the  time  of  arrival  data  to  define  the  other  two 
layers  as  before.  Calculation  PD16  also  assumed  that  the 
boundary  between  the  working  point  layer  and  the  layer  above 
it  was  the  location  of  the  water  table  as  well.  The  effective 
stress  law  was  only  used  below  the  water  table. 

Calculation  PD16  has  been  the  most  encouraging  calcula¬ 
tion  made  to  this  date,  giving  good  agreement  with  much  of  the 
data.  However,  the  coarse  zoning  in  this  preliminary  calcula¬ 
tion  has  not  allowed  us  to  model  with  sufficient  accuracy  the 
spall  at  the  free  surface.  Future  calculations  will  locate 
the  water  table  more  accurately  and  attempt  to  model  the  spall 
in  greater  detail. 

A  parallel  effort,  using  one-dimensional  calculations, 
is  underway  to  model  the  French  nuclear  tests  in  the  granite 
of  the  Hoggar  Massif  in  the  Sahara  Desert.  Laboratory  tests 
on  core  samples  have  indicated  that  the  strength  and  other 
material  properties  of  this  rock  are  fairly  similar  to  NTS 
granodiorite.  Of  course,  the  Hoggar  rock  is  very  dry.  To 
first  approximation  (material  properties  are  not  identical) , 
calculations  in  Hoggar  granite  may  then  be  made  with  the  same 
material  properties  as  for  PILEDRIVER,  but  excluding  the  effec¬ 
tive  stress  law.  Analogous  calculations  to  PILEDRIVER  Calcula¬ 
tions  1  and  2  have  been  made  for  Hoggar  without  the  effective 
stress  model.  These  calculations  show  considerably  smaller 
cube  root  scaled  cavity  radii  than  PILEDRIVER,  but  are  within 
three  percent  and  nine  percent,  respectively,  of  the  cavity 
radii  reported  by  the  French.  The  RVP  transforms  computed 


31 


srsrrMS.  socncc  an o  torrwAmc 


from  the  Hoggar  calculations  show  relatively  little  peaking  of 
the  spectra.  The  static  RDP  is  approximately  a  factor  of  three 
to  five  lower  than  computed  for  PILEDRIVER. 

In  summary,  partial  saturation  coupled  with  an  effective 
stress  law  provides  a  constitutive  model  which  gives  an 
acceptable  fit  to  the  PILEDRIVER  free  field  ground  motion 
data  and  the  measured  cavity  radius.  A  match  to  the  reported 
Hoggar  cavity  radii  is  obtained  by  assuming  strengh  character¬ 
istics  appropriate  for  a  dry  rock  environment.  The  resulting 
seismic  coupling  is  a  factor  of  three  to  five  greater  for  the 
PILEDRIVER  environment  than  for  Hoggar.  If  we  believe  the 
reported  Hoggar  cavities,  then  for  identical  propagation  path 
characteristics,  a  given  yield  in  the  PILEDRIVER  environment 
would  generate  teleseismic  body  waves  which  would  be  three 
to  five  times  larger  than  those  produced  by  the  same  yield 
in  the  Hoggar  environment. 
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VI.  DECOUPLING  CALCULATIONS  IN  SALT 


Figure  8  shows  the  amplitude  spectra  of  the  reduced 
velocity  potential  (RVP)  for  SALMON  free  field  data  (Springer, 
et  al. ,  1968).  We  have  calculated  RVP  spectra  for  0.38  kt 
(the  STERLING  yield)  detonated  in  a  18  meter  cavity. 

Two  calculations  were  completed  for  this  cavity  radius, 
one  involving  detonation  in  a  mined  cavity  and  the  other 
detonation  in  a  cavity  with  the  rock  outside  the  cavity 
weakened  from  shock  heating.  The  amount  of  thermal  weakening 
was  estimated  from  a  calculation  of  the  SALMON  event. 

Figure  9  shows  the  decoupling  ratios  between  the  RVP 
spectra  shown  in  Figure  8  (scaled  to  0.38  kt)  and  these  two 
calculations.  For  frequencies  less  than  7  Hz  the  decoupling 
ratio  is  approximately  130  for  the  mined  cavity  and  110  for 
the  thermally  weakened  SALMON  cavity.  Figure  10  shows  the 
decoupling  ratio  between  the  RVP  spectrum  in  Figure  8  and  the 
spectrum  from  5.3  kt  detonated  in  a  50  meter  mined  cavity. 

The  decoupling  ratio  is  approximately  130  but  the  corner  fre¬ 
quency  has  decreased  to  2.5  Hz  due  to  detonation  in  the  larger 
cavi ty . 

Near  Source  Attenuation 

Figure  11  shows  a  comparison  between  calculated  and 
observed  peak  radial  stress  and  particle  velocity  from  the 
SALMON  event.  The  elastic  radius  from  the  calculation  occurs 
at  approximately  300  meters  while  the  data  shows  no  apparent 
change  of  slope  over  the  entire  distance  range. 

A  possible  explanation  for  this  discrepancy  is  that 
the  peak  particle  velocity  (and  stress)  is  sensitive  to  the 
high  frequency  components  of  the  seismic  source  while  the 
nonlinear  constitutive  models  for  fracture,  plastic  flow. 
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FREQUENCY 

Figure  9.  Decoupling  ratio  between  SALMON  d< 
calculations  for  0.38  kt  detonatec 


3  CL 


irreversible  pore  collapse  and  effective  stress  are  not  in¬ 
tended  to  model  the  anelastic  attenuation  of  small  deformation 
elastic  waves. 

In  order  to  quantify  this  hypothesis  we  have  investi¬ 
gated  the  effect  of  Q  on  peak  displacement, particle  velocity 
and  radial  stress. 

Figure  12  shows  the  change  in  the  calculated  RVP 
spectrum  for  SALMON  over  the  distance  range  0.3km<R<0.6km 
due  to  convolution  with  the  attenuation  operator. 

irf  (R  -  0.3  km) 
aQ 

where  Q  *  10.  The  change  in  the  spectrum  for  frequencies  less 
than  10  Hz  is  not  dramatic  even  with  attenuation  this  severe. 

Figure  13  shows  the  effect  of  Q  on  peak  displacement, 
particle  velocity  and  stress.  We  see  a  large  change  in  the 
attenuation  rate  for  peak  velocity  and  stress  when  Q  changes 
from  infinity  to  ten.  This  is  a  lesser  change  in  the  attenua¬ 
tion  of  peak  displacement,  almost  certainly  caused  by  the 
lower  frequencies  responsible  for  this  quantity. 

It  seems  that  most  of  the  problems  associated  with 
the  observed  attenuation  rate  of  peak  velocity,  stress  and 
displacement  versus  those  predicted  by  nonlinear  finite 
difference  codes  can  be  resolved  by  the  addition  of  anelastic 
attenuation  to  the  equivalent  source  predicted  by  the  codes. 
This  has  been  our  approach  from  the  beginning.  In  addition, 
measurements  of  the  peaks  can  be  used  to  infer  the  Q  ap¬ 
propriate  for  propagation  of  this  source  to  distances  beyond 
the  elastic  radius. 
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Figure  12.  The  effect  of  Q  -  10  on  the  RVP  spectra. 
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Figure  13.  The  effect  of  Q  on  peak  displacement,  stress  and 
velocity. 
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